# ------------------ importing simulation library ---------------------- #
from actinNetwork import *
from numpy import arccos,sqrt,median,argmax,std,unique,savetxt,vstack
from numpy.random import normal,choice
from scipy.special import erf
from scipy.stats import zscore
from matplotlib.pyplot import subplots,bar,xticks,axhline,setp
import csv
# display plots in notebook
%matplotlib inline
# ------------------------------- parameters for simulation ----------------------------- #
# initial number of filaments N, and monomer size dN, and range in nm
N = 200; dN = 5.5; filRange = 1000.0;
# running time T, precision dt, and frame rate ds
T = 10.0; dt = 0.001; ds = 2.0;
# velocity and width of polymerisation profile
Fo = 0.2; dw = 5*dN
vo = 250.0
# ------------------------------- velocity pertubation ---------------------------------- #
# time and duration of perturbation
tau = 4.0; dtau = 0.1
# velocity change during perturbation
dF = 0.2
# velocity perturbation as heaviside function centered at tau
def F(t) :
return Fo - dF * heavisidePi( t-tau-sqrt(dtau), 2.0*sqrt(dtau) )
# ------------------------------- rate assumptions ---------------------------------- #
# define polymerisation rate in terms of
# critical angle given unperturbed velocity
phiCrit = 50.0 * pi/180
# rates for travelling heavisides
lambdaRate = vo/(dN*cos(phiCrit))
betaRate = 0.5
kappaRate = 1000.0
# defining rate functions
def rLambda(x,y,theta,t) :
return lambdaRate*heavisidePi(y,dw)
def rBeta(x,y,theta,t) :
return betaRate*heavisidePi(y,dw)
def rKappa(x,y,theta,t) :
return kappaRate*heavisideTheta(-y-dw/2)+1.0*heavisidePi(y,dw)
plot([t for t in arange(0,T,dt)],[F(t) for t in arange(0,T,dt)])
plot([t for t in arange(0,T,dt)],[Fo + exp(-(t-1.065*tau)**2) for t in arange(0,T,dt)])
# binning cross x direction
xEdges1 = linspace(0,2000,27)
xCenters1 = (xEdges1[:-1] + xEdges1[1:]) / 2.0
xEdges2 = linspace(0,2000,14)
xCenters2 = (xEdges2[:-1] + xEdges2[1:]) / 2.0
# binning across phi direction
phiEdges1 = array([0,20,50,90])*pi/180
phiCenters1 = ( phiEdges1[1:]+phiEdges1[:-1] ) / 2.0
phiEdges2 = linspace(-pi/2,pi/2,10)
phiCenters2 = ( phiEdges2[1:]+phiEdges2[:-1] ) / 2.0
# function that calculates relevant statistics given a network object that has already been run
def getStatistics(networkObject,xEdges1,xEdges2,phiEdges1,phiEdges2) :
# getting positions
xFil = networkObject.getPositions( networkObject.Monomers )
xBranch = networkObject.getPositions( networkObject.Branches )
xCap = networkObject.getPositions( networkObject.Caps )
# angular orientation of filaments
phiFil = networkObject.getAngles ( networkObject.Monomers )
phiBranch = networkObject.getAngles ( networkObject.Branches )
phiCap = networkObject.getAngles ( networkObject.Caps )
# calculating normalised filament denisty in x direction
filamentDensity = histogram( xFil.T[1], bins = xEdges1)[0]
filamentDensity = filamentDensity/float(max(filamentDensity))
# calculating zscore branching and capping rates in x direction
branchRate = zscore( histogram( xBranch.T[1], bins = xEdges1)[0] )
capRate = zscore( histogram( xCap.T[1], bins = xEdges1)[0] )
# binning monomers into x-slices
filMasks = [ ( x <= xFil.T[1] ) * ( xFil.T[1] <= y ) for x,y in zip(xEdges2[:-1],xEdges2[1:]) ]
# and using these bins to calculate angule histograms at each slice
filNumber = array([ histogram( abs(unique(phiFil[mask])), bins = phiEdges1)[0] for mask in filMasks ]).astype(float)
phiDensity = array([ histogram( phiFil[mask], bins=phiEdges2)[0] for mask in filMasks ])
# calculation of order parameter
orderEdges = array([-50,-30,-10,10,30,50])*pi/180
orderHist = [ histogram(phiFil[mask],bins=orderEdges)[0] for mask in filMasks ]
orderParam = [ ( N[2]-(N[0]+N[-1])/2.0 ) / float( N[2]+(N[0]+N[-1])/2.0 ) for N in orderHist ]
# return relevant statistics
return filamentDensity,branchRate,capRate,filNumber,phiDensity,orderParam
#xSteady = vstack([n.getPositions(n.Frontier).T[0],uniform(low=0.0, high=dw/2.0, size=len(n.Frontier))]).T;
#dxSteady = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
# ------------------------------- initialisation ---------------------------------- #
# initial distribution of network
try : # using previous positions and angles if possible
# getting previous positions and angles
#xInit = n.getPositions(n.Frontier) - array(len(n.Frontier)*[[0.0,p(n.tElapsed)]])
#dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
xInit = xSteady
dxInit = dxSteady
except : # otherwise initialise uniformly
# sampling from uniform distribution given N
thetaDist = uniform(-pi,pi,size=N)
xDist = uniform(low=0.0, high=filRange, size=N)
yDist = uniform(low=0.0, high=dw/2.0, size=N)
# passing samples to inital variables
xInit = array([ [x,y] for x,y in zip(xDist,yDist) ])
dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in thetaDist ])
# initialise network object
n = network(rLambda, rBeta, rKappa,
xSeed = xInit, dxSeed = dxInit,
branchSigma = 15.0*pi/180.0,
forceDirection = True, recordHistory = True)
# ------------------------------- simulation ---------------------------------- #
# evolve network up to specified time
#n.exportData(dt,ds,n.tElapsed+T,Fext=1.0e-1)
#n.evolve(dt,n.tElapsed+T,Fext=1.0)
# while networks grows and up until time tFinal
while n.nBarbed != 0 and n.tElapsed <= T :
# evolve
n.timeStep(dt, Fext = F(n.tElapsed) )
### ------------------------------- Filament Plot ------------------------------ ###
# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )
figure(figsize=(20,5))
# plotting all monomers as different colour points
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=4)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=4)
# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18)
ylabel(r"Distance, $y$ / nm",fontsize=18)
ylim(0,n.xBoundary); xlim(0,2000);
# ------------------------------- calculate statistics ---------------------------------- #
filamentDensity,branchRate,capRate,filNumber,phiDist,orderParam = getStatistics(n,xEdges1,xEdges2,phiEdges1,phiEdges2)
### ------------------------------- Density Plot ------------------------------ ###
f,ax1 = subplots()
ax2 = ax1.twinx()
f.set_figheight(5)
f.set_figwidth(20)
# plotting rates and density
ax1.plot( xCenters1, filamentDensity, '#00ff00',mec='#00ff00',marker="o",mew="4",ms=12,linewidth=4)
ax2.plot( xCenters1, branchRate, '#2737ff',mec='#2737ff',marker="o",mew="4",ms=12,linewidth=4)
ax2.plot( xCenters1, capRate, '#ff0000',mec='#ff0000',marker="o",mew="4",ms=12,linewidth=4)
# plotting options
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
ax1.set_ylabel(r"Percentage Filament Denisty, $\rho(x)$",fontsize=16)
ax2.set_ylabel(r"Standardised Rate, $\omega(x)$", fontsize=16)
ax1.set_ylim(0,1); ax1.set_xlim(0,)
### ------------------------------- Capping Ratios ------------------------------ ###
figure(figsize=(20,5))
# plot capping ratios
plot( xCenters2, filNumber.T[0],'k',
mec='k',marker="o",mew="4",ms=12,linewidth=4,
label=r"0$^{\circ}$- 20$^{\circ}$")
plot( xCenters2, 2.0/3.0*filNumber.T[1],'orange',
mec='orange',marker="o",mew="4",ms=12,linewidth=4,
label=r"20$^{\circ}$- 50$^{\circ}$")
plot( xCenters2, 1.0/2.0*filNumber.T[2],'0.7',
mec='0.7', marker="o",mew="4",ms=12,linewidth=4,
label=r"50$^{\circ}$- 90$^{\circ}$")
# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=16); l = legend(title=r"Angle $\phi$",fontsize=16)
ylabel(r"Filament Number $N(x\vert\phi)$",fontsize=16)
setp(l.get_title(),fontsize=16)
xlim(0,)
### ------------------------------- Angle Order ------------------------------ ###
f,ax1 = subplots()
ax2 = ax1.twinx()
f.set_figheight(5)
f.set_figwidth(20)
# plotting angle distributions
for x,y in zip([ linspace(x-40,x+40,len(phiCenters2)) for x in xCenters2 ],phiDist) :
ax1.bar( x, y, width=10, color="gray", edgecolor="k", linewidth=2, align="center")
# and order parameter
ax2.plot( xCenters2, orderParam, "r.-", marker="o", mec="r",mew="4",ms=12,linewidth=4, label="Order Parameter")
ax2.axhline( 0, color="k", linewidth=1);
# plotting options
ax1.set_ylabel(r"Filament Density, $\rho(\phi\vert x)$",fontsize=16);
ax2.set_ylabel(r"Order Parameter,$\alpha(x)$",fontsize=16);
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
xlim(0,);
# initialising arrays to gather results
FD = []
BR = []
CR = []
AD = []
OP = []
# number of times to run the simulation
nRuns = 20
# looping over nRuns
for run in range(nRuns) :
# ------------------------------- initialisation ---------------------------------- #
# initial distribution of network
try : # using previous positions and angles if possible
# getting previous positions and angles
# xInit = n.getPositions(n.Frontier) - array(len(n.Frontier)*[[0.0,p(n.tElapsed)]])
# dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
#xSteady = vstack([n.getPositions(n.Frontier).T[0],uniform(low=0.0, high=dw/2.0, size=len(n.Frontier))]).T;
#dxSteady = array([ [dN*sin(theta),dN*cos(theta)] for theta in n.getAngles(n.Frontier) ])
xInit = xSteady
dxInit = dxSteady
except : # otherwise initialise uniformly
# sampling from uniform distribution given N
thetaDist = uniform(-pi,pi,size=N)
xDist = uniform(low=0.0, high=filRange, size=N)
yDist = uniform(low=0.0, high=dw/2.0, size=N)
# passing samples to inital variables
xInit = array([ [x,y] for x,y in zip(xDist,yDist) ])
dxInit = array([ [dN*sin(theta),dN*cos(theta)] for theta in thetaDist ])
# initialise network object
n = network(rLambda, rBeta, rKappa,
xSeed = xInit, dxSeed = dxInit,
branchSigma = 15.0*pi/180.0,
forceDirection = True, recordHistory = True)
# ------------------------------- simulation ---------------------------------- #
# while networks grows and up until time tFinal
while n.nBarbed != 0 and n.tElapsed <= T :
# evolve
n.timeStep(dt, Fext = F(n.tElapsed) )
### ----------------------------- Processing Data --------------------------- ###
filamentDensity,branchRate,capRate,filNumber,phiDist,orderParam = getStatistics(n,xEdges1,xEdges2,phiEdges1,phiEdges2)
FD += [filamentDensity]
BR += [branchRate]
CR += [capRate]
AD += [filNumber]
OP += [orderParam]
### ------------------------------- Filament Plot ------------------------------ ###
# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )
figure(figsize=(20,5))
# plotting all monomers as different colour points
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=5,alpha=0.2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=6,alpha=0.5)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=6,alpha=0.5)
# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18)
ylabel(r"Distance, $y$ / nm",fontsize=18)
ylim(0,n.xBoundary); xlim(0,2000);
### ------------------------------- Density Plot ------------------------------ ###
f,ax1 = subplots()
ax2 = ax1.twinx()
f.set_figheight(5)
f.set_figwidth(20)
# plotting rates and density
ax1.errorbar( xCenters1, mean(array(FD),axis=0),yerr=std(array(FD),axis=0), color='#00ff00',mec='#00ff00',marker="o",mew="4",ms=0,linewidth=4)
ax2.errorbar( xCenters1, mean(array(BR),axis=0),yerr=std(array(BR),axis=0), color='#2737ff',mec='#2737ff',marker="o",mew="4",ms=0,linewidth=4)
ax2.errorbar( xCenters1, mean(array(CR),axis=0),yerr=std(array(CR),axis=0), color='#ff0000',mec='#ff0000',marker="o",mew="4",ms=0,linewidth=4)
# plotting options
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=18)
ax1.set_ylabel(r"Percentage Filament Denisty, $\rho(x)$",fontsize=18)
ax2.set_ylabel(r"Standardised Rate, $\omega(x)$", fontsize=16)
ax1.set_ylim(0,1); ax2.set_ylim(-5,5); ax1.set_xlim(0,2000);
f,ax = subplots()
f.set_figheight(5)
f.set_figwidth(20)
# plot capping ratios
ax.errorbar( xCenters2, mean(array(AD),axis=0).T[0],yerr=std(array(AD),axis=0).T[0],color='k',
mec='k',marker="o",mew="4",ms=0,linewidth=4,
label=r"0$^{\circ}$- 20$^{\circ}$")
ax.errorbar( xCenters2, 2.0/3.0*mean(array(AD),axis=0).T[1],yerr=2.0/3.0*std(array(AD),axis=0).T[1],color='orange',
mec='orange',marker="o",mew="4",ms=0,linewidth=4,
label=r"20$^{\circ}$- 50$^{\circ}$")
ax.errorbar( xCenters2, 1.0/2.0*mean(array(AD),axis=0).T[2],yerr=1.0/2.0*std(array(AD),axis=0).T[2],color='0.7',
mec='0.7', marker="o",mew="4",ms=0,linewidth=4,
label=r"50$^{\circ}$- 90$^{\circ}$")
# plotting options
xlabel(r"Distance, $x$ / nm",fontsize=18); legend(fontsize=16)
ylabel(r"Percentage Angular Denisty $\rho(x\vert\phi)$",fontsize=18)
xlim(0,2000);
f,ax1 = subplots()
ax2 = ax1.twinx()
f.set_figheight(5)
f.set_figwidth(20)
# plotting angle distributions
for x,y in zip([ linspace(x-80,x+80,len(phiCenters2)) for x in xCenters2 ],phiDist) :
ax1.bar( x, y, width=20, color="none", edgecolor="k", linewidth=2, align="center")
# and order parameter
ax2.errorbar( xCenters2, mean(array(OP),axis=0),yerr=std(array(OP),axis=0), color="r", linewidth=4, label="Order Parameter")
ax2.axhline( 0, color="k", linewidth=1);
# plotting options
ax1.set_ylabel(r"Number of Monomers, $N(\phi\vert x)$",fontsize=16);
ax2.set_ylabel(r"Order Parameter, $\alpha(x)$",fontsize=16);
ax1.set_xlabel(r"Distance, $x$ / nm",fontsize=16)
xlim(0,2000);
# ---------------------------------------save everything!! ----------------------- #
savetxt("X-FD-BR-CR-correlated.csv",
array([xCenters1,
mean(array(FD),axis=0),std(array(FD),axis=0),
mean(array(BR),axis=0),std(array(BR),axis=0),
mean(array(CR),axis=0),std(array(CR),axis=0) ]).T,delimiter=",")
savetxt("X-AD-correlated.csv",array([xCenters2,
1.0/1.0*mean(array(AD),axis=0).T[0],1.0/1.0*std(array(AD),axis=0).T[0],
2.0/3.0*mean(array(AD),axis=0).T[1],2.0/3.0*std(array(AD),axis=0).T[1],
1.0/2.0*mean(array(AD),axis=0).T[2],1.0/2.0*std(array(AD),axis=0).T[2] ]).T,delimiter=",")
savetxt("X-OP-correlated.csv",array([xCenters2,
mean(array(OP),axis=0),std(array(OP),axis=0) ]).T,delimiter=",")
# binning monomers into x-slices
phiFil = n.getAngles ( n.Monomers )
xFil = n.getPositions( n.Monomers )
filMasks = [ ( x <= xFil.T[1] ) * ( xFil.T[1] <= y ) for x,y in zip(xEdges2[:-1],xEdges2[1:]) ]
b = open("SAD-correlated.csv", 'w')
a = csv.writer(b)
data = [[str(x*180.0/pi) for x in phiFil[mask]] for mask in filMasks ]
a.writerows(data)
b.close()
for i,x in enumerate(OP[:2]) :
plot(x,color=str(i/float(len(OP))))
axhline( 0, color="k", linewidth=1)
### ------------------------------- Filament Plot ------------------------------ ###
# getting positions
xFil = n.getPositions( n.Monomers )
xBranch = n.getPositions( n.Branches )
xCap = n.getPositions( n.Caps )
figure(figsize=(80,20))
# plotting all monomers as different colour points
plot(xFil.T[1],xFil.T[0],'g',marker=".",linewidth=0,ms=12,alpha=0.2)
plot(xBranch.T[1],xBranch.T[0],'#2737ff',marker=".",linewidth=0,ms=15,alpha=0.5)
plot(xCap.T[1],xCap.T[0],'#ff0000',marker=".",linewidth=0,ms=15,alpha=0.5)
# plotting options
axis("off");
n.nBarbed